library(tidyverse)
library(janitor)
library(viridis)
library(fst)
library(lfe)
library(broom)
library(stargazer)
library(glue)
library(nnet)
library(fastDummies)

options(stringsAsFactors = F,
        dplyr.summarise.inform = F)

# BEGIN ANALYSIS ----------------------------------------------------------

vf <- read_fst("vf_for_analysis_anon.fst")

turnout_vars = c("general_2018", "primary_2018",
                 "general_2016", "primary_2016")

turnout_labs = list("general_2018"="2018 General",
                    "primary_2018"="2018 Primary",
                    "general_2016"="2016 General",
                    "primary_2016"="2016 Primary")

auto_vars = c("auto_in_hh", "drivers_license")

auto_labs = list("auto_in_hh" = "Auto in HH",
                 "drivers_license" = "Drivers License")

covar_order_list <- c("auto_in_hh", "drivers_license", "gender", "white", "age")

small_felm <- function(m){
  m$fe <- NULL
  #m$residuals <- NULL
  #m$r.residuals <- NULL
  m$fitted.values <- NULL
  m$response <- NULL
  m$cfactor <- NULL
  m$inv <- NULL
  m$STATS$promotion <- NULL
  m$clustervar <- NULL
  m
}

# Raw Differences ---------------------------------------------------------

for(v in turnout_vars) {
  print(v)
  writeLines(as.character(round(mean(pull(vf, !!v),na.rm=T)*100,0)),
             glue("tables/mean_turnout_{v}.tex"),sep = "%")
  
  auto_summ <- vf %>% 
    group_by(auto_in_hh) %>%
    summarize(x=mean(!!sym(v),na.rm=T))
  writeLines(as.character(round(auto_summ$x[auto_summ$auto_in_hh==0]*100,0)),
             glue("tables/mean_turnout_nocar_{v}.tex"),sep = "%")
  writeLines(as.character(round(auto_summ$x[auto_summ$auto_in_hh==1]*100,0)),
             glue("tables/mean_turnout_car_{v}.tex"),sep = "%")
  writeLines(as.character(round((auto_summ$x[auto_summ$auto_in_hh==1]-
                                   auto_summ$x[auto_summ$auto_in_hh==0])*100,0)),
             glue("tables/diff_turnout_carvnocar_{v}.tex"),sep = "%")
  
  auto_summ_plot <- ggplot(auto_summ) + 
    geom_histogram(aes(x=auto_in_hh,y=x),stat="identity") + 
    geom_text(aes(x=auto_in_hh,y=x+0.03,label=paste0(round(100*x,0),"%"))) + 
    scale_x_continuous("Auto in Household",breaks=c(0,1),labels=c("No","Yes")) + 
    scale_y_continuous(glue("{turnout_labs[[v]]} Turnout"),labels=scales::percent_format(),limits=c(0,1)) + 
    theme_minimal()
  ggsave(auto_summ_plot, filename = glue("figures/auto_summary_{v}.pdf"),
         height=4,width=4, units="in", device=cairo_pdf)
}

x <- map_dfr(turnout_vars, function(v) {
  auto_summ <- vf %>% 
    group_by(auto_in_hh) %>%
    summarize(x=mean(!!sym(v),na.rm=T)) %>%
    mutate(lab=turnout_labs[[v]])
})

x %>% filter(str_detect(lab, "2018")) %>%
  ggplot() +
  geom_histogram(aes(x=auto_in_hh,y=x),stat="identity") + 
  geom_text(aes(x=auto_in_hh,y=x+0.03,
                label=paste0(round(100*x,0),"%")),
            family="Arial Narrow", size=3) + 
  scale_x_continuous("Auto in Household",breaks=c(0,1),labels=c("No","Yes")) + 
  facet_wrap(~lab, ncol=2) +
  scale_y_continuous("Turnout",labels=scales::percent_format(),limits=c(0,1)) + 
  theme_minimal(base_family = "Arial Narrow") +
  theme(strip.background = element_rect(fill="grey"))
ggsave(filename = glue("figures/auto_summary_2018b.pdf"),
       height=3,width=4, units="in", device=cairo_pdf)
ggsave(filename = glue("figures/auto_summary_2018b.png"),
       height=3,width=4, units="in")

# Fixed Effects Models ----------------------------------------------------

for(a in auto_vars) {
  m <- list()
  for(v in turnout_vars) {
    for(fe in c("0", "county", "county + precinct", "address_id")) {
      print(glue("{v} {fe} {a}"))
      #f = as.formula(paste0(v, " ~  ", a, " | ", fe, " | 0 | 0"))
      #m[[glue("{v}_{fe}_{a}")]] <- felm(f, data=vf) %>% small_felm()
      f = as.formula(paste0(v, " ~  ", a, " + gender + white + age | ", fe, " | 0 | 0"))
      m[[glue("{v}_{fe}_{a}_cov")]] <- felm(f, data=vf) %>% small_felm()
    }
  }
  for(y in c("2016", "2018")) {
    stargazer(m[str_detect(names(m), y) & str_detect(names(m), a)],
              keep.stat = c("n","rsq","adj.rsq"),
              dep.var.labels = c(glue("{y} General Turnout"),glue("{y} Primary Turnout")),
              column.labels = NULL,
              covariate.labels = c(
                auto_labs[[a]],
                "Male",
                "White",
                "Age"
              ),
              order=covar_order_list,
              add.lines = list(c("FE for County",rep(c("", "$\\checkmark$", "", ""),2, each=1)),
                               c("FE for Precinct",rep(c("", "", "$\\checkmark$", ""),2, each=1)),
                               c("FE for Address",rep(c("", "", "", "$\\checkmark$"),2, each=1))),
              # notes.append = T, notes="Standard errors clustered by county",
              out = glue("tables/turnout_fe_{y}_{a}.tex"),
              float = F,
              star.cutoffs = c(.01, NA, NA),
              notes = "$^{*}$p$<$0.01", notes.append=F
    )
  }
}

remove(m)

# Travel Time -------------------------------------------------------------

travel_sample <- filter(vf, voter_in_randomsample==1) %>%
  mutate(traveltime_diff_car = seconds_min_noncar - seconds_min_car,
         traveltime_min_diff_car = traveltime_diff_car/60)

ggplot(travel_sample) + 
  geom_density(aes(x=minutes_min_car),fill="black",col="black",alpha=0.8) + 
  geom_density(aes(x=minutes_min_noncar),fill="gray50",col="black",alpha=0.6) + 
  annotate(geom = "text",col="black",x=8,y=0.15,label="With Car", hjust=0, family="Arial Narrow", size=3) + 
  annotate(geom = "text",col="black",x=45,y=0.015,label="Without Car", hjust=0, family="Arial Narrow", size=3) + 
  scale_x_continuous("Travel time to polling place (min.)",breaks=c(15,45,75,105,135),limits=c(0,160)) + 
  scale_y_continuous("Density") + 
  # ggtitle("All MI Counties") + 
  theme_minimal(base_family = "Arial Narrow")
ggsave(filename = "figures/traveltime_carv nocar_randomsample_density.pdf",height=3,width=5, device = cairo_pdf)

writeLines(as.character(round(mean(travel_sample$traveltime_min_diff_car,na.rm=T),1)),"tables/avg_traveltimediff_mins.tex",sep = "%")
writeLines(as.character(round(median(travel_sample$traveltime_min_diff_car,na.rm=T),1)),"tables/median_traveltimediff_mins.tex",sep = "%")

ggplot(travel_sample) + 
  geom_density(aes(x=traveltime_min_diff_car),fill="darkgray",col="black",alpha=1) + 
  scale_x_continuous("Difference in travel time to polling place without car vs. with car (min.)",breaks=c(15,45,75,105,135),limits=c(-5,160)) + 
  scale_y_continuous("Density") + 
  # ggtitle("All MI Counties") + 
  theme_minimal(base_family = "Arial Narrow")
ggsave(filename = "figures/traveltimediff_carvnocar_randomsample_density.pdf",height=3,width=5, device = cairo_pdf)

## Descriptive diffs by travel time:
n_groups = 4
travel_sample <- travel_sample %>%
  mutate(traveltime_nocar_bin = ntile(minutes_min_noncar,n = n_groups))
x1 <- travel_sample %>% 
  group_by(quartile=traveltime_nocar_bin, auto_in_hh) %>% 
  summarize(x=mean(general_2018, na.rm=T)) %>%
  drop_na()

x <- quantile(travel_sample$minutes_min_noncar,probs = c(0,0.25,0.5,0.75,1),na.rm = T) %>% 
  round(digits=1)
quartile_labs <- data.frame(quartile = c(1:4),
                            labs = c(glue("1st Quartile\n(<{x[2]} min.)"),
                                     glue("2nd Quartile\n({x[2]}-{x[3]} min.)"),
                                     glue("3rd Quartile\n({x[3]}-{x[4]} min.)"),
                                     glue("4th Quartile\n(>{x[4]} min.)")))
x1$labs <- quartile_labs$labs[match(x1$quartile,quartile_labs$quartile)]

avg <- mean(travel_sample$general_2018, na.rm=T)

ggplot(x1, aes(x=x, y=reorder(labs, desc(labs)))) +
  geom_vline(xintercept = avg,lty=2,col="red") + 
  geom_line(aes(group=labs)) +
  geom_point(aes(fill=as.factor(auto_in_hh)), pch=21) +
  geom_text(data=filter(x1, auto_in_hh==0),
            aes(x=x-0.03,label=paste0(round(100*x,0), "%")),hjust="right", family="Arial Narrow", size=3) + 
  geom_text(data=filter(x1, auto_in_hh==1),
            aes(x=x+0.03,label=paste0(round(100*x,0), "%")),hjust="left", family="Arial Narrow", size=3) + 
  annotate("text",x=avg-0.01,y=2.5,
           label=paste0("1% sample average\nturnout = ",round(avg,2)*100,"%"),hjust="right",col="red",size=3, family="Arial Narrow") +
  scale_fill_manual(values=c("white", "black"),
                    labels=c("No Auto in HH", "Auto in HH"), name="") +
  labs(y="", x="Turnout in 2018 general election") +
  scale_x_continuous(labels=scales::percent_format(), limits=c(0, 1)) +
  theme_minimal(base_family = "Arial Narrow") +
  theme(strip.background = element_rect(fill="grey"),
        legend.position = "bottom",
        legend.background = element_rect(color="black", size=.25),
        legend.margin = margin(1,10,1,1))
ggsave(filename = "figures/average_turnout_traveltimequartiles.pdf",width=4,height=3, device=cairo_pdf)

# Travel Time: Quartile Models --------------------------------------------

tt_overall <- felm(general_2018 ~  auto_in_hh + gender + white + age | county + precinct | 0 | 0, 
                  data=travel_sample)
n_groups = 4
travel_sample <- travel_sample %>%
  mutate(traveltime_diff_car_bin = ntile(traveltime_min_diff_car,n = n_groups))

m <- list()
for(i in 1:n_groups) {
  print(i)
  m[[i]] <- felm(general_2018 ~  auto_in_hh + gender + white + age | county + precinct | 0 | 0, 
                 data=filter(travel_sample, traveltime_diff_car_bin==i))
  
}

stargazer(m,
          keep.stat = c("n","rsq","adj.rsq"),
          dep.var.labels = "2018 Turnout",
          column.labels = c("1st Quartile","2nd Quartile","3rd Quartile","4th Quartile"),
          covariate.labels = c(
            "Auto in HH",
            # "\\# autos in HH",
            "Male",
            "White",
            "Age"
          ),
          add.lines = list(c("FE for Precinct",rep("$\\checkmark$", n_groups))),
          # notes.append = T, notes="Standard errors clustered by county",
          out = "tables/general_2018_travel_time_quantiles.tex",
          float = F,
          star.cutoffs = c(.01, NA, NA),
          notes = "$^{*}$p$<$0.01", notes.append=F
)

coefs <- map(m, ~ tidy(., conf.int = T, conf.level = .95)[1,]) %>% bind_rows(.id="q")
x <- quantile(travel_sample$traveltime_min_diff_car,probs = c(0,0.25,0.5,0.75,1),na.rm = T) %>% round(digits=1)
coefs$subset <- c(glue("1st Quartile\n(<{x[2]} min.)"),
                  glue("2nd Quartile\n({x[2]}-{x[3]} min.)"),
                  glue("3rd Quartile\n({x[3]}-{x[4]} min.)"),
                  glue("4th Quartile\n(>{x[4]} min.)"))

ggplot(coefs, aes(y=reorder(subset, desc(q)), x=estimate, xmin=conf.low, xmax=conf.high)) +
  geom_vline(xintercept = 0,lty=2) +
  geom_vline(xintercept = (tt_overall$coefficients[1]),lty=2,col="red") + 
  geom_errorbarh(height=.1) +
  geom_point() +
  geom_text(aes(x=conf.low-0.005,label=round(estimate,3)),hjust="right", family="Arial Narrow", size=3) + 
  annotate("text",x=tt_overall$coefficients[1]-0.005,y=2.5,
           label=paste0("1% sample average\neffect = ",round(tt_overall$coefficients[1],3)),hjust="right",col="red",size=3, family="Arial Narrow") +
  labs(y="Difference in travel time to polls\nwithout car-with a car",
       x="Within-precinct effect of auto in household on\nprob. turning out in 2018 general election (95% CIs)") + 
  theme_minimal(base_family = "Arial Narrow")
ggsave(filename = "figures/travel_time_quartiles_general_2018.pdf",width=4,height=3, device=cairo_pdf)
writeLines(as.character(round(coefs$estimate[1]*100,1)),"tables/travel_time_quartiles_general_2018_q1.tex",sep = "%")
writeLines(as.character(round(coefs$estimate[4]*100,1)),"tables/travel_time_quartiles_general_2018_q4.tex",sep = "%")

rm(m)



# Travel Times: Quartiles for Matched Addresses ---------------------------

d <- filter(vf, !is.na(address_id))
tt_overall <- felm(general_2018 ~  auto_in_hh + gender + white + age | address_id | 0 | 0, 
                   data=d)
n_groups = 4
d <- d %>% 
  mutate(traveltime_diff_car = seconds_min_noncar - seconds_min_car,
                  traveltime_min_diff_car = traveltime_diff_car/60) %>%
  mutate(traveltime_diff_car_bin = ntile(traveltime_min_diff_car,n = n_groups))
#quantile(travel_sample$traveltime_diff_car,probs = c(0,0.25,0.5,0.75,1),na.rm = T)

m <- list()
for(i in 1:n_groups) {
  print(i)
  m[[i]] <- felm(general_2018 ~  auto_in_hh + gender + white + age | address_id | 0 | 0, 
                 data=filter(d, traveltime_diff_car_bin==i))
  
}

stargazer(m,
          keep.stat = c("n","rsq","adj.rsq"),
          dep.var.labels = "2018 Turnout",
          column.labels = c("1st Quartile","2nd Quartile","3rd Quartile","4th Quartile"),
          covariate.labels = c(
            "Auto in HH",
            # "\\# autos in HH",
            "Male",
            "White",
            "Age"
          ),
          add.lines = list(c("FE for Address",rep("$\\checkmark$", n_groups))),
          # notes.append = T, notes="Standard errors clustered by county",
          out = "tables/general_2018_travel_time_quantiles_address.tex",
          float = F,
          star.cutoffs = c(.01, NA, NA),
          notes = "$^{*}$p$<$0.01", notes.append=F
)

coefs <- map(m, ~ tidy(., conf.int = T, conf.level = .95)[1,]) %>% bind_rows(.id="q")
x <- quantile(travel_sample$traveltime_min_diff_car,probs = c(0,0.25,0.5,0.75,1),na.rm = T) %>% round(digits=1)
coefs$subset <- c(glue("1st Quartile\n(<{x[2]} min.)"),
                  glue("2nd Quartile\n({x[2]}-{x[3]} min.)"),
                  glue("3rd Quartile\n({x[3]}-{x[4]} min.)"),
                  glue("4th Quartile\n(>{x[4]} min.)"))

ggplot(coefs, aes(y=reorder(subset, desc(q)), x=estimate, xmin=conf.low, xmax=conf.high)) +
  geom_vline(xintercept = 0,lty=2) +
  geom_vline(xintercept = (tt_overall$coefficients[1]),lty=2,col="red") + 
  geom_errorbarh(height=.1) +
  geom_point() +
  geom_text(aes(x=conf.low-0.005,label=round(estimate,3)),hjust="right", family="Arial Narrow", size=3) + 
  annotate("text",x=tt_overall$coefficients[1]-0.005,y=2.5,
           label=paste0("1% sample average\neffect = ",round(tt_overall$coefficients[1],3)),hjust="right",col="red",size=3, family="Arial Narrow") +
  labs(y="Difference in travel time to polls\nwithout car-with a car",
       x="Within-address effect of auto in household on\nprob. turning out in 2018 general election (95% CIs)") + 
  theme_minimal(base_family = "Arial Narrow")
ggsave(filename = "figures/travel_time_quartiles_general_2018_address.pdf",width=4,height=3, device=cairo_pdf)
writeLines(as.character(round(coefs$estimate[1]*100,1)),"tables/travel_time_quartiles_general_2018_q1_address.tex",sep = "%")
writeLines(as.character(round(coefs$estimate[4]*100,1)),"tables/travel_time_quartiles_general_2018_q4_address.tex",sep = "%")


# Race & Age Effects ------------------------------------------------------
## Car access descriptives
auto_access_race_summ <- vf %>% 
  group_by(race) %>%
  summarize(v=mean(auto_in_hh,na.rm=T))

writeLines(as.character(round(auto_access_race_summ$v[auto_access_race_summ$race=="White"]*100,0)),"tables/mean_car_white.tex",sep = "%")
writeLines(as.character(round(auto_access_race_summ$v[auto_access_race_summ$race=="Black"]*100,0)),"tables/mean_car_black.tex",sep = "%")
writeLines(as.character(round(auto_access_race_summ$v[auto_access_race_summ$race=="Hispanic"]*100,0)),"tables/mean_car_hispanic.tex",sep = "%")

(auto_access_race_summ_plot <- ggplot(auto_access_race_summ) + 
    geom_histogram(aes(x=race,y=v),stat="identity") + 
    geom_text(aes(x=race,y=v+0.03,label=paste0(round(100*v,0),"%"))) + 
    scale_x_discrete("Race/ethnicity") + 
    scale_y_continuous("Auto in household",labels=scales::percent_format(),limits=c(0,1)) + 
    theme_minimal())
ggsave(auto_access_race_summ_plot,filename = "figures/auto_access_byrace.pdf",height=4,width=4)

auto_access_age_summ <- vf %>% 
  filter(!is.na(age_cat)) %>%
  group_by(age_cat) %>%
  summarize(v=mean(auto_in_hh,na.rm=T)) %>%
  mutate(age_cat = str_replace(age_cat," and over","\nand over"))

writeLines(as.character(round(auto_access_age_summ$v[auto_access_age_summ$age_cat=="18-24"]*100,0)),"tables/mean_car_1824.tex",sep = "%")
writeLines(as.character(round(auto_access_age_summ$v[auto_access_age_summ$age_cat=="25-34"]*100,0)),"tables/mean_car_2534.tex",sep = "%")
writeLines(as.character(round(auto_access_age_summ$v[auto_access_age_summ$age_cat=="45-54"]*100,0)),"tables/mean_car_4554.tex",sep = "%")
writeLines(as.character(round(auto_access_age_summ$v[auto_access_age_summ$age_cat=="55-64"]*100,0)),"tables/mean_car_5564.tex",sep = "%")
writeLines(as.character(round(auto_access_age_summ$v[auto_access_age_summ$age_cat=="65 years\nand over"]*100,0)),"tables/mean_car_65plus.tex",sep = "%")

(auto_access_age_summ_plot <- ggplot(auto_access_age_summ) + 
    geom_histogram(aes(x=age_cat,y=v),stat="identity") + 
    geom_text(aes(x=age_cat,y=v+0.03,label=paste0(round(100*v,0),"%"))) + 
    scale_x_discrete("Age") + 
    scale_y_continuous("Auto in household",labels=scales::percent_format(),limits=c(0,1)) + 
    theme_minimal())
ggsave(auto_access_age_summ_plot,filename = "figures/auto_access_byage.pdf",height=4,width=4)


m <- list()
for(i in c("White", "Black", "Hispanic", "Asian", "Other")) {
  print(i)
  m[[paste0("race_", i)]] <- felm(general_2018 ~ auto_in_hh + 
                            age + gender | county + precinct | 0 | 0, 
                            data=filter(vf,race==i))
}
for(i in c("18-24", "25-34", "35-44", "45-54", "55-64", "65 years and over")){
  print(i)
  m[[paste0("age_", i)]] <- felm(general_2018 ~ auto_in_hh + 
                                    white + gender | county + precinct | 0 | 0, 
                                  data=filter(vf,age_cat==i))
}

stargazer(m[str_detect(names(m), "race")],
          keep.stat = c("n","rsq","adj.rsq"),
          dep.var.labels = c(glue("2018 General Turnout")),
          column.labels = c("White", "Black", "Hispanic", "Asian", "Other"),
          covariate.labels = c(
            "Auto in HH",
            "Male",
            "Age"
          ),
          order=covar_order_list,
          add.lines = list(c("FE for Precinct",rep("$\\checkmark$", 5))),
          # notes.append = T, notes="Standard errors clustered by county",
          out = glue("tables/turnout_by_subgroup_race.tex"),
          float = F,
          star.cutoffs = c(.01, NA, NA),
          notes = "$^{*}$p$<$0.01", notes.append=F
)

stargazer(m[str_detect(names(m), "age")],
          keep.stat = c("n","rsq","adj.rsq"),
          dep.var.labels = c(glue("2018 General Turnout")),
          column.labels = c("18-24", "25-34", "35-44", "45-54", "55-64", "65+"),
          covariate.labels = c(
            "Auto in HH",
            "Male",
            "White"
          ),
          order=covar_order_list,
          add.lines = list(c("FE for Precinct",rep("$\\checkmark$", 6))),
          # notes.append = T, notes="Standard errors clustered by county",
          out = glue("tables/turnout_by_subgroup_age.tex"),
          float = F,
          star.cutoffs = c(.01, NA, NA),
          notes = "$^{*}$p$<$0.01", notes.append=F
)

coefs <- map(m, ~ tidy(., conf.int = T, conf.level = .95)[1,]) %>% 
  bind_rows(.id="q") %>%
  separate(q, into=c("cat", "lab"), sep="_") %>%
  mutate(lab=ifelse(lab=="65 years and over", "65+", lab),
         cat=str_to_sentence(cat))

ggplot(coefs, aes(y=factor(lab, levels = rev(lab)), x=estimate, xmin=conf.low, xmax=conf.high)) +
  geom_vline(xintercept = 0,lty=2) +
  geom_errorbarh(height=.1) +
  geom_point() +
  geom_text(aes(x=conf.low-0.005,label=round(estimate,3)),hjust="right", family="Arial Narrow", size=3) + 
  facet_wrap(~cat, ncol = 1, scale="free_y") +
  labs(y="",
       x="Within-precinct effect of auto in household on\nprob. turning out in 2018 general election (95% CIs)") + 
  theme_minimal(base_family = "Arial Narrow") +
  theme(strip.background = element_rect(fill="grey"))
ggsave(filename = "figures/fe_auto_in_hh_age_race.pdf",width=4,height=4, device=cairo_pdf)


x1 <- vf %>% group_by(lab=race, auto_in_hh) %>% summarize(x=mean(general_2018, na.rm=T)) %>%
  mutate(cat="Race")
x2 <- vf %>% group_by(lab=age_cat, auto_in_hh) %>% summarize(x=mean(general_2018, na.rm=T)) %>%
  mutate(lab=ifelse(lab=="65 years and over", "65+", lab)) %>%
  drop_na() %>%
  mutate(cat="Age")
x <- bind_rows(x1, x2)
avg <- mean(vf$general_2018, na.rm=T)

  ggplot(x, aes(x=x, y=reorder(lab, desc(lab)))) +
  geom_line(aes(group=lab)) +
  geom_point(aes(fill=as.factor(auto_in_hh)), pch=21) +
  geom_text(data=filter(x, auto_in_hh==0),
            aes(x=x-0.03,label=paste0(round(100*x,1), "%")),hjust="right", family="Arial Narrow", size=3) + 
    geom_text(data=filter(x, auto_in_hh==1),
              aes(x=x+0.03,label=paste0(round(100*x,1), "%")),hjust="left", family="Arial Narrow", size=3) + 
  scale_fill_manual(values=c("white", "black"),
                    labels=c("No Auto in HH", "Auto in HH"), name="") +
  facet_wrap(~cat, ncol = 1, scale="free_y") +
    labs(y="", x="Turnout in 2018 general election") +
    scale_x_continuous(labels=scales::percent_format(), limits=c(0, 1)) +
  theme_minimal(base_family = "Arial Narrow") +
  theme(strip.background = element_rect(fill="grey"),
        legend.position = "bottom",
        legend.background = element_rect(color="black", size=.25),
        legend.margin = margin(1,10,1,1))
ggsave(filename = "figures/average_turnout_age_race.pdf",width=4,height=4, device=cairo_pdf)

y <- x %>% filter(cat=="Race") %>%
  pivot_wider(names_from=auto_in_hh, values_from=x, names_prefix="auto")
for(i in 1:nrow(y)) {
  writeLines(as.character(round(y$auto1[i]*100,1)),
             glue("tables/mean_turnout_car_{str_to_lower(y$lab[i])}.tex"),
             sep = "%")
  writeLines(as.character(round(y$auto0[i]*100,1)),
             glue("tables/mean_turnout_nocar_{str_to_lower(y$lab[i])}.tex"),
             sep = "%")
}

writeLines(as.character(round((y$auto0[1]-y$auto0[2])*100,1)),
           glue("tables/diff_turnout_nocar_white_black.tex"),
           sep = "%")
writeLines(as.character(round((y$auto1[1]-y$auto1[2])*100,1)),
           glue("tables/diff_turnout_car_white_black.tex"),
           sep = "%")

writeLines(as.character(round((y$auto0[1]-y$auto0[3])*100,1)),
           glue("tables/diff_turnout_nocar_white_hispanic.tex"),
           sep = "%")
writeLines(as.character(round((y$auto1[1]-y$auto1[3])*100,1)),
           glue("tables/diff_turnout_car_white_hispanic.tex"),
           sep = "%")


# Vote Method -------------------------------------------------------------

mode_general_2018 <- vf %>% 
  filter(!is.na(address_id)) %>%
  filter(!is.na(state_gen_2018_mode)) %>%
  group_by(auto_in_hh,state_gen_2018_mode) %>%
  summarize(n=n()) %>%
  group_by(auto_in_hh) %>%
  mutate(total = sum(n),
         perc = n/total) %>%
  mutate(lab="2018 General", m=state_gen_2018_mode)
writeLines(as.character(round(mode_general_2018$perc[mode_general_2018$auto_in_hh==0 & mode_general_2018$state_gen_2018_mode=="ABS"]*100,0)),"tables/rate_abs_state18_matchedHH_nocar.tex",sep = "%")
writeLines(as.character(round(mode_general_2018$perc[mode_general_2018$auto_in_hh==1 & mode_general_2018$state_gen_2018_mode=="ABS"]*100,0)),"tables/rate_abs_state18_matchedHH_car.tex",sep = "%")
writeLines(as.character(round(mode_general_2018$perc[mode_general_2018$auto_in_hh==0 & mode_general_2018$state_gen_2018_mode=="ELEC"]*100,0)),"tables/rate_inpers_state18_matchedHH_nocar.tex",sep = "%")
writeLines(as.character(round(mode_general_2018$perc[mode_general_2018$auto_in_hh==1 & mode_general_2018$state_gen_2018_mode=="ELEC"]*100,0)),"tables/rate_inpers_state18_matchedHH_car.tex",sep = "%")
writeLines(as.character(round((mode_general_2018$perc[mode_general_2018$auto_in_hh==1 & mode_general_2018$state_gen_2018_mode=="ELEC"]-mode_general_2018$perc[mode_general_2018$auto_in_hh==0 & mode_general_2018$state_gen_2018_mode=="ELEC"])*100,0)),"tables/diff_inpers_state18_matchedHH_carvnocar.tex",sep = "%")
writeLines(as.character(round((mode_general_2018$perc[mode_general_2018$auto_in_hh==1 & mode_general_2018$state_gen_2018_mode=="ELEC"]-mode_general_2018$perc[mode_general_2018$auto_in_hh==0 & mode_general_2018$state_gen_2018_mode=="ELEC"])/mode_general_2018$perc[mode_general_2018$auto_in_hh==0 & mode_general_2018$state_gen_2018_mode=="ELEC"]*100,0)),"tables/diff_perc_inpers_state18_matchedHH_carvnocar.tex",sep = "%")

mode_primary_2018 <- vf %>% 
  filter(!is.na(address_id)) %>%
  filter(!is.na(state_pri_2018_mode)) %>%
  group_by(auto_in_hh,state_pri_2018_mode) %>%
  summarize(n=n()) %>%
  group_by(auto_in_hh) %>%
  mutate(total = sum(n),
         perc = n/total) %>%
  mutate(lab="2018 Primary", m=state_pri_2018_mode)

bind_rows(mode_general_2018, mode_primary_2018) %>%
ggplot() + 
    geom_histogram(aes(x=auto_in_hh,y=perc,fill=m),stat="identity",position=position_dodge(0.9)) + 
    geom_text(aes(x=auto_in_hh,y=perc+0.05,label=paste0(round(100*perc,0),"%"),
                  group=m),position=position_dodge(0.9),
              family="Arial Narrow", size=3) + 
  facet_wrap(~lab, ncol=1, scales="free_x") +
    scale_x_continuous(breaks=c(0,1),labels=c("No Auto in HH","Auto in HH")) + 
    scale_y_continuous(labels=scales::percent_format(),limits=c(0,1)) + 
  labs(x="", y="") +
    scale_fill_viridis_d("",aesthetics=c("fill","color"),limits=c("ABS","ELEC","NO"),labels=c("Absentee","In-person","Did not vote")) + 
    theme_minimal(base_family = "Arial Narrow") +
  theme(legend.position = c(0.8,0.9)) +
  theme(strip.background = element_rect(fill="grey"))
ggsave(filename = "figures/vote_mode_2018.pdf",height=6,width=4, device=cairo_pdf)

last_plot() + facet_wrap(~lab, ncol=2, scales="free_x") +
  theme(legend.position = c(0.4,0.85))
ggsave(filename = "figures/vote_mode_2018_horiz.pdf",height=3,width=6, device=cairo_pdf)
ggsave(filename = "figures/vote_mode_2018_horiz.png",height=3,width=6)

x <- filter(vf, !is.na(address_id)) %>%
  mutate(abs_gen_2018=as.numeric(state_gen_2018_mode=="ABS"),
         abs_pri_2018=as.numeric(state_pri_2018_mode=="ABS"),
         elec_gen_2018=as.numeric(state_gen_2018_mode=="ELEC"),
         elec_pri_2018=as.numeric(state_pri_2018_mode=="ELEC"))

m <- list()
m[["abs_gen_2018"]] <- felm(abs_gen_2018 ~ auto_in_hh + gender + white + age | address_id | 0 | 0,
     data=x)
m[["elec_gen_2018"]] <- felm(elec_gen_2018 ~ auto_in_hh + gender + white + age | address_id | 0 | 0,
                            data=x)
m[["abs_pri_2018"]] <- felm(abs_pri_2018 ~ auto_in_hh + gender + white + age | address_id | 0 | 0,
                            data=x)
m[["elec_pri_2018"]] <- felm(elec_pri_2018 ~ auto_in_hh + gender + white + age | address_id | 0 | 0,
                            data=x)

stargazer(m, keep.stat = c("n","rsq","adj.rsq"),
          dep.var.labels = c("2018 General Absentee",
                             "2018 General In-Person",
                             "2018 Primary Absentee",
                             "2018 Primary In-Person"),
          column.labels = NULL,
          covariate.labels = c(
            "Auto in HH",
            "Male",
            "White",
            "Age"
          ),
          order=covar_order_list,
          add.lines = list(c("FE for Address",rep("$\\checkmark$", 4))),
          # notes.append = T, notes="Standard errors clustered by county",
          out = glue("tables/vote_mode_within_address_general_2018.tex"),
          float = F,
          star.cutoffs = c(.01, NA, NA),
          notes = "$^{*}$p$<$0.01", notes.append=F
)

remove(m)

## Vote Method - MNL --------------------------------------------------------
vf <- vf %>%
  mutate(gen_2018_choice = relevel(factor(state_gen_2018_mode),ref = "NO"))

mode_fit <- (multinom(gen_2018_choice ~ auto_in_hh + gender + white + age,
                       data=vf))
summary(mode_fit)
stargazer(mode_fit,
          type="text",
          out = "tables/choice_2018_mnl.txt",
          dep.var.labels = c("Choose absentee over not voting", "Choosing in-person over not voting"),
          #df = F,omit.stat = c("rsq","ser"),
          covariate.labels = c("Car in household","Male","White","Age")
)

# Predicted probabilities:
vf$male <- as.numeric(vf$gender=="M")
vf <- vf %>%
  fastDummies::dummy_cols("age_cat") %>%
  rename(age_1824 = `age_cat_18-24`,
         age_2534 = `age_cat_25-34`,
         age_3544 = `age_cat_35-44`,
         age_4554 = `age_cat_45-54`,
         age_5564 = `age_cat_55-64`,
         age_65plus = `age_cat_65 years and over`)

mode_fit_dummies <- (multinom(gen_2018_choice ~ auto_in_hh + male + white + age_1824 + age_2534 + age_3544 + age_4554 + age_5564 + age_65plus,
                       data=filter(vf,!is.na(address_id) & !is.na(state_gen_2018_mode))))
summary(mode_fit_dummies)

D.car <- model.frame(mode_fit_dummies)
D.car$auto_in_hh <- 1
pred.car <- predict(mode_fit_dummies, newdata = D.car, type = "probs")
car.probs <- apply(pred.car,2,mean)
D.nocar <- model.frame(mode_fit_dummies)
D.nocar$auto_in_hh <- 0
pred.nocar <- predict(mode_fit_dummies, newdata = D.nocar, type = "probs")
nocar.probs <- apply(pred.nocar,2,mean)

D.male <- model.frame(mode_fit_dummies)
D.male$male <- 1
pred.male <- predict(mode_fit_dummies, newdata = D.male, type = "probs")
male.probs <- apply(pred.male,2,mean)
D.female <- model.frame(mode_fit_dummies)
D.female$male <- 0
pred.female <- predict(mode_fit_dummies, newdata = D.female, type = "probs")
female.probs <- apply(pred.female,2,mean)

D.white <- model.frame(mode_fit_dummies)
D.white$white <- 1
pred.white <- predict(mode_fit_dummies, newdata = D.white, type = "probs")
white.probs <- apply(pred.white,2,mean)
D.nonwhite <- model.frame(mode_fit_dummies)
D.nonwhite$white <- 0
pred.nonwhite <- predict(mode_fit_dummies, newdata = D.nonwhite, type = "probs")
nonwhite.probs <- apply(pred.nonwhite,2,mean)

D.1824 <- model.frame(mode_fit_dummies)
D.1824$age_1824 <- 1
D.1824$age_2534 <- 0
D.1824$age_3544 <- 0
D.1824$age_4554 <- 0
D.1824$age_5564 <- 0
D.1824$age_65plus <- 0
pred.1824 <- predict(mode_fit_dummies, newdata = D.1824, type = "probs")
age1824.probs <- apply(pred.1824,2,mean)
D.2534 <- model.frame(mode_fit_dummies)
D.2534$age_1824 <- 0
D.2534$age_2534 <- 1
D.2534$age_3544 <- 0
D.2534$age_4554 <- 0
D.2534$age_5564 <- 0
D.2534$age_65plus <- 0
pred.2534 <- predict(mode_fit_dummies, newdata = D.2534, type = "probs")
age2534.probs <- apply(pred.2534,2,mean)
D.3544 <- model.frame(mode_fit_dummies)
D.3544$age_1824 <- 0
D.3544$age_2534 <- 0
D.3544$age_3544 <- 1
D.3544$age_4554 <- 0
D.3544$age_5564 <- 0
D.3544$age_65plus <- 0
pred.3544 <- predict(mode_fit_dummies, newdata = D.3544, type = "probs")
age3544.probs <- apply(pred.3544,2,mean)
D.4554 <- model.frame(mode_fit_dummies)
D.4554$age_1824 <- 0
D.4554$age_2534 <- 0
D.4554$age_3544 <- 0
D.4554$age_4554 <- 1
D.4554$age_5564 <- 0
D.4554$age_65plus <- 0
pred.4554 <- predict(mode_fit_dummies, newdata = D.4554, type = "probs")
age4554.probs <- apply(pred.4554,2,mean)
D.5564 <- model.frame(mode_fit_dummies)
D.5564$age_1824 <- 0
D.5564$age_2534 <- 0
D.5564$age_3544 <- 0
D.5564$age_4554 <- 0
D.5564$age_5564 <- 1
D.5564$age_65plus <- 0
pred.5564 <- predict(mode_fit_dummies, newdata = D.5564, type = "probs")
age5564.probs <- apply(pred.5564,2,mean)
D.65plus <- model.frame(mode_fit_dummies)
D.65plus$age_1824 <- 0
D.65plus$age_2534 <- 0
D.65plus$age_3544 <- 0
D.65plus$age_4554 <- 0
D.65plus$age_5564 <- 0
D.65plus$age_65plus <- 1
pred.65plus <- predict(mode_fit_dummies, newdata = D.65plus, type = "probs")
age65plus.probs <- apply(pred.65plus,2,mean)



## Predicted Probabilities ##
prob_tab <- data.frame(cbind(c("Auto in HH","No Auto in HH","Male","Female","White","Non-white","Age: 18-24","Age: 25-34","Age: 35-44","Age: 45-54","Age: 55-64","Age: 65+"),
                             rbind(car.probs,nocar.probs,male.probs,female.probs,white.probs,nonwhite.probs,age1824.probs,age2534.probs,age3544.probs,age4554.probs,age5564.probs,age65plus.probs))) %>%
  mutate(NO = round(100*as.numeric(NO),2),ABS = round(100*as.numeric(ABS),2),ELEC = round(100*as.numeric(ELEC),2)) %>%
  rename(Variable = V1, `Did not vote` = NO, Absentee = ABS, `In-person` = ELEC)

  
write_csv(prob_tab,
          file = "tables/mode_2018_mnl_pred_probs.txt")

print(xtable::xtable(prob_tab),include.rownames=F,file="tables/mode_2018_mnl_pred_probs.tex",floating=F,booktabs=T)

## Coefficients:
stargazer(mode_fit_dummies,
          out = "tables/mode_2018_mnl_dummies.tex",
          dep.var.labels = c("Choose absentee over not voting", "Choose in-person over not voting"),
          #df = F,omit.stat = c("rsq","ser"),
          covariate.labels = c("Auto in HH","Male","White","Age: 18-24","Age: 25-34","Age: 35-44","Age: 45-54","Age: 55-64","Age: 65$+$"),
          float = F
)

# County Results ----------------------------------------------------------

for(v in turnout_vars) {
  print(v)
  cty_res <- vf %>% select(county, v=!!v, auto_in_hh, gender, white, age) %>% 
    nest(data=-county) %>%
    mutate(m = map(data, ~ lm(v ~ auto_in_hh + gender + white + age, 
                              data = .)),
           results = map(m, ~ tidy(., conf.int=T))) %>% 
    unnest(results) %>% 
    select(-data, -m) %>% 
    filter(term=="auto_in_hh")
  
  cty_res %>%
    ggplot(aes(x=reorder(str_to_sentence(county), -estimate), y=estimate, ymin=conf.low, ymax=conf.high)) +
    geom_errorbar() +
    geom_point() +
    scale_y_continuous(limits=c(0, .5)) +
    labs(x="", y="Effect of Car Ownership on Turnout") +
    theme_minimal(base_family = "Arial Narrow") +
    theme(axis.text.x = element_text(angle=45, hjust=1, size=5))
  ggsave(glue("figures/county_regressions_{v}.pdf"), width=6, height=3, 
         units="in", device = cairo_pdf)
}


#  Auto x License Interaction ---------------------------------------------
## Licenses - descriptives
license_race_summ <- vf %>% 
  group_by(race) %>%
  summarize(v=mean(drivers_license,na.rm=T))

writeLines(as.character(round(license_race_summ$v[license_race_summ$race=="White"]*100,0)),"tables/mean_license_white.tex",sep = "%")
writeLines(as.character(round(license_race_summ$v[license_race_summ$race=="Black"]*100,0)),"tables/mean_license_black.tex",sep = "%")
writeLines(as.character(round(license_race_summ$v[license_race_summ$race=="Hispanic"]*100,0)),"tables/mean_license_hispanic.tex",sep = "%")

(license_race_summ_plot <- ggplot(license_race_summ) + 
    geom_histogram(aes(x=race,y=v),stat="identity") + 
    geom_text(aes(x=race,y=v+0.03,label=paste0(round(100*v,0),"%"))) + 
    scale_x_discrete("Race/ethnicity") + 
    scale_y_continuous("Matched to drivers'\nlicense data",labels=scales::percent_format(),limits=c(0,1.03)) + 
    theme_minimal())
ggsave(license_race_summ_plot,filename = "figures/license_byrace.pdf",height=4,width=4)

license_age_summ <- vf %>% 
  filter(!is.na(age_cat)) %>%
  group_by(age_cat) %>%
  summarize(v=mean(drivers_license,na.rm=T)) %>%
  mutate(age_cat = str_replace(age_cat," and over","\nand over"))

writeLines(as.character(round(license_age_summ$v[license_age_summ$age_cat=="18-24"]*100,0)),"tables/mean_license_1824.tex",sep = "%")
writeLines(as.character(round(license_age_summ$v[license_age_summ$age_cat=="25-34"]*100,0)),"tables/mean_license_2534.tex",sep = "%")
writeLines(as.character(round(license_age_summ$v[license_age_summ$age_cat=="45-54"]*100,0)),"tables/mean_license_4554.tex",sep = "%")
writeLines(as.character(round(license_age_summ$v[license_age_summ$age_cat=="55-64"]*100,0)),"tables/mean_license_5564.tex",sep = "%")
writeLines(as.character(round(license_age_summ$v[license_age_summ$age_cat=="65 years\nand over"]*100,0)),"tables/mean_license_65plus.tex",sep = "%")

(license_age_summ_plot <- ggplot(license_age_summ) + 
    geom_histogram(aes(x=age_cat,y=v),stat="identity") + 
    geom_text(aes(x=age_cat,y=v+0.03,label=paste0(round(100*v,0),"%"))) + 
    scale_x_discrete("Age") + 
    scale_y_continuous("Matched to drivers'\nlicense data",labels=scales::percent_format(),limits=c(0,1.03)) + 
    theme_minimal())
ggsave(license_age_summ_plot,filename = "figures/license_byage.pdf",height=4,width=4)


m <- list()
for(v in turnout_vars) {
  for(fe in c("0", "county", "county + precinct", "address_id")) {
    print(glue("{v} {fe}"))
    #f = as.formula(paste0(v, " ~  ", a, " | ", fe, " | 0 | 0"))
    #m[[glue("{v}_{fe}_{a}")]] <- felm(f, data=vf) %>% small_felm()
    f = as.formula(paste0(v, " ~ auto_in_hh*drivers_license + gender + white + age | ", fe, " | 0 | 0"))
    m[[glue("{v}_{fe}_cov")]] <- felm(f, data=vf) %>% small_felm()
  }
}
for(y in c("2016", "2018")) {
  stargazer(m[str_detect(names(m), y)],
            keep.stat = c("n","rsq","adj.rsq"),
            dep.var.labels = c(glue("{y} General Turnout"),glue("{y} Primary Turnout")),
            column.labels = NULL,
            covariate.labels = c(
              "Auto in HH",
              "Drivers License",
              "Auto in HH x Drivers License",
              "Male",
              "White",
              "Age"
            ),
            order=c("auto_in_hh", "drivers_license", "auto_in_hh:drivers_license",
                    "gender", "white", "age"),
            add.lines = list(c("FE for County",rep(c("", "$\\checkmark$", "", ""),2, each=1)),
                             c("FE for Precinct",rep(c("", "", "$\\checkmark$", ""),2, each=1)),
                             c("FE for Address",rep(c("", "", "", "$\\checkmark$"),2, each=1))),
            # notes.append = T, notes="Standard errors clustered by county",
            out = glue("tables/turnout_fe_{y}_auto_x_drivers.tex"),
            float = F,
            star.cutoffs = c(.01, NA, NA),
            notes = "$^{*}$p$<$0.01", notes.append=F
  )
}

remove(m)

# L2 Covariates Models ----------------------------------------------------

for(a in auto_vars) {
  m <- list()
  for(v in turnout_vars) {
    for(fe in c("0", "county", "county + precinct", "address_id")) {
      print(glue("{v} {fe} {a}"))
      #f = as.formula(paste0(v, " ~  ", a, " | ", fe, " | 0 | 0"))
      #m[[glue("{v}_{fe}_{a}")]] <- felm(f, data=vf) %>% small_felm()
      f = as.formula(paste0(v, " ~  ", a, " + gender + white + age + 
                            l2_est_hh_income +
                            l2_education + l2_home_owner_or_renter | ", fe, " | 0 | 0"))
      m[[glue("{v}_{fe}_{a}_cov")]] <- felm(f, data=vf) %>% small_felm()
    }
  }
  for(y in c("2016", "2018")) {
    stargazer(m[str_detect(names(m), y) & str_detect(names(m), a)],
              keep.stat = c("n","rsq","adj.rsq"),
              dep.var.labels = c(glue("{y} General Turnout"),glue("{y} Primary Turnout")),
              column.labels = NULL,
              covariate.labels = c(
                auto_labs[[a]],
                "Male",
                "White",
                "Age",
                "Est. HH Income",
                "HS Diploma",
                "Vocational Degree",
                "Some College",
                "College Degree",
                "Grad Degree",
                "Renter"
              ),
              order=c("auto_in_hh", "drivers_license",
                      "gender", "white", "age",
                      "l2_est_hh_income",
                      "l2_education",
                      "l2_home_owner_or_renter"),
              add.lines = list(c("FE for County",rep(c("", "$\\checkmark$", "", ""),2, each=1)),
                               c("FE for Precinct",rep(c("", "", "$\\checkmark$", ""),2, each=1)),
                               c("FE for Address",rep(c("", "", "", "$\\checkmark$"),2, each=1))),
              # notes.append = T, notes="Standard errors clustered by county",
              out = glue("tables/turnout_fe_{y}_{a}_l2_covars.tex"),
              float = F,
              star.cutoffs = c(.01, NA, NA),
              notes = "$^{*}$p$<$0.01", notes.append=F
    )
  }
}

remove(m)


# 2020 Models -------------------------------------------------------------

m <- list()
for(a in auto_vars) {
  v <- "general_2020"
    for(fe in c("0", "county2020", "county2020 + precinct2020")) {
      print(glue("{v} {fe} {a}"))
      #f = as.formula(paste0(v, " ~  ", a, " | ", fe, " | 0 | 0"))
      #m[[glue("{v}_{fe}_{a}")]] <- felm(f, data=vf) %>% small_felm()
      f = as.formula(paste0(v, " ~  ", a, " + gender + white + age | ", fe, " | 0 | 0"))
      m[[glue("{v}_{fe}_{a}_cov")]] <- felm(f, data=vf) %>% small_felm()
    }
}

y <- "2020"
stargazer(m,
          keep.stat = c("n","rsq","adj.rsq"),
          dep.var.labels = c(glue("{y} General Turnout")),
          column.labels = NULL,
          covariate.labels = c(
            auto_labs[[1]],
            auto_labs[[2]],
            "Male",
            "White",
            "Age"
          ),
          order=covar_order_list,
          add.lines = list(c("FE for County",rep(c("", "$\\checkmark$", ""),2, each=1)),
                           c("FE for Precinct",rep(c("", "", "$\\checkmark$"),2, each=1))),
          # notes.append = T, notes="Standard errors clustered by county",
          out = glue("tables/turnout_fe_{y}.tex"),
          float = F,
          star.cutoffs = c(.01, NA, NA),
          notes = "$^{*}$p$<$0.01", notes.append=F
)

remove(m)


# Voters with Drivers Licenses Only ---------------------------------------

dl <- filter(vf, drivers_license==1)
a = "auto_in_hh"

  m <- list()
  for(v in turnout_vars) {
    for(fe in c("0", "county", "county + precinct", "address_id")) {
      print(glue("{v} {fe} {a}"))
      #f = as.formula(paste0(v, " ~  ", a, " | ", fe, " | 0 | 0"))
      #m[[glue("{v}_{fe}_{a}")]] <- felm(f, data=vf) %>% small_felm()
      f = as.formula(paste0(v, " ~  ", a, " + gender + white + age | ", fe, " | 0 | 0"))
      m[[glue("{v}_{fe}_{a}_cov")]] <- felm(f, data=dl) %>% small_felm()
    }
  }
  for(y in c("2016", "2018")) {
    stargazer(m[str_detect(names(m), y) & str_detect(names(m), a)],
              keep.stat = c("n","rsq","adj.rsq"),
              dep.var.labels = c(glue("{y} General Turnout"),glue("{y} Primary Turnout")),
              column.labels = NULL,
              covariate.labels = c(
                auto_labs[[a]],
                "Male",
                "White",
                "Age"
              ),
              order=covar_order_list,
              add.lines = list(c("FE for County",rep(c("", "$\\checkmark$", "", ""),2, each=1)),
                               c("FE for Precinct",rep(c("", "", "$\\checkmark$", ""),2, each=1)),
                               c("FE for Address",rep(c("", "", "", "$\\checkmark$"),2, each=1))),
              # notes.append = T, notes="Standard errors clustered by county",
              out = glue("tables/turnout_fe_dl_only_{y}_{a}.tex"),
              float = F,
              star.cutoffs = c(.01, NA, NA),
              notes = "$^{*}$p$<$0.01", notes.append=F
    )
  }

remove(m)


# Survey Data Evidence ---------------------------------------
## TAPS data collected by Wash U
## data available here: https://wc.wustl.edu/taps-data-archive
## below are merged files with profile data
dat14 <- read_rds("taps_survey_2014.rds")
dat15 <- read_rds("taps_survey_2015.rds")

## Data cleaning and recodes ---------------------------------------
dat14 <- dat14 %>%
  mutate(
    owncar14 = dplyr::recode(car16s37,"yes"=1,"no"=0,"Refused"=NA_real_),
    vote2014 = recode(sprres7s36,"yes" = 1, "no" = 0, "Refused" = NA_real_)
  )

dat15 <- dat15 %>%
  mutate(
    owncar15 = dplyr::recode(car16s42,"yes"=1,"no"=0,"Refused"=NA_real_),
    vote2014 = recode(sprres7s47,"yes" = 1, "no" = 0, "Refused" = NA_real_)
  )

## Analyses  ---------------------------------------
fit_control_14 <- lm(vote2014 ~ owncar14 + 
                       factor(ppethm) + 
                       factor(ppgender) + 
                       factor(ppeducat) + 
                       factor(ppagect4),
                     data=dat14)

fit_control_15 <- lm(vote2014 ~ owncar15 + 
                       factor(ppethm) + 
                       factor(ppgender) + 
                       factor(ppeducat) + 
                       factor(ppagect4),
                     data=dat15)

stargazer(fit_control_14,fit_control_15,
          dep.var.labels = "Reported Voting in Nov. 2014",
          column.labels = c("Nov. 2014 Survey","Oct. 2015 Survey"),
          covariate.labels = c("Reported driving a car regularly, Dec. 2014",
                               "Reported driving a car regularly, May 2015",
                               "Race/Eth.: Black, non-Hispanic",
                               "Race/Eth.: Other, non-Hispanic",
                               "Race/Eth.: Hispanic",
                               "Race/Eth.: 2+ Races, non-Hispanic",
                               "Female",
                               "Education: High school degree",
                               "Education: Some college",
                               "Education: Bachelor's degree or higher",
                               "Age: 30-44",
                               "Age: 45-59",
                               "Age: 60+"),
          omit.stat = c("ser","adj.rsq"),
          float=F,
          notes=c("Omitted category for race is White, non-Hispanic",
                  "Omitted category for education is Less than high school",
                  "Omitted category for age is 18-29"),
          out = "tables/turnout_survey_taps.tex"
)

